module glc_override_frac

#include "shr_assert.h"

  !---------------------------------------------------------------------------
  ! !DESCRIPTION:
  ! This module provides functionality to allow overriding the ice fractions (i.e., the
  ! ice_covered field) and topographic heights that are sent to the coupler

  ! !USES:

  use glc_kinds_mod
  
  implicit none
  private
  save

  ! !PUBLIC MEMBER FUNCTIONS:
  public :: init_glc_frac_overrides ! initialize stuff in this module, including reading the namelist
  public :: frac_overrides_enabled  ! return true if overrides are enabled, false otherwise
  public :: do_frac_overrides       ! do all overrides

  ! !PRIVATE MEMBER FUNCTIONS:
  private :: read_namelist        ! read namelist setting options in this module
  private :: apply_increase_frac
  private :: apply_decrease_frac
  private :: apply_rearrange_freq
  private :: time_since_baseline  ! return time (days) since the baseline for adjustments
  
  
  ! !PRIVATE DATA MEMBERS:
  logical(log_kind) :: enable_frac_overrides     ! whether the overrides in this module are enabled for this run
  integer(int_kind) :: decrease_override_delay   ! time delay before beginning decrease_frac overrides (days)
  integer(int_kind) :: increase_override_delay   ! time delay before beginning increase_frac overrides (days)
  integer(int_kind) :: rearrange_override_delay  ! time delay before beginning rearrange_freq overrides (days)
  real(r8)          :: decrease_frac             ! fractional decrease per day (should be positive)
  real(r8)          :: increase_frac             ! fractional increase per day
  integer(int_kind) :: rearrange_freq            ! frequency (days) at which we rearrange elevation classes

  ! Assumed maximum topographic height. It's okay for heights to go above this value, but
  ! they may not be handled exactly as desired by the overrides here.
  real(r8), parameter :: max_height = 3500._r8
  !---------------------------------------------------------------------------

contains

  !-----------------------------------------------------------------------
  subroutine init_glc_frac_overrides
    !
    ! !DESCRIPTION:
    ! Initialize stuff in this module, including reading the namelist
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    !
    ! !LOCAL VARIABLES:

    character(len=*), parameter :: subname = 'init_glc_frac_overrides'
    !-----------------------------------------------------------------------
    
    call read_namelist

  end subroutine init_glc_frac_overrides
  

  !-----------------------------------------------------------------------
  subroutine read_namelist
    !
    ! !DESCRIPTION:
    ! Read namelist setting options in this module
    !
    ! !USES:
    use glc_files      , only : nml_filename
    use glc_constants  , only : stdout, nml_in, blank_fmt, ndelim_fmt
    use glc_communicate, only : my_task, master_task
    use glc_broadcast  , only : broadcast_scalar
    use glc_exit_mod   , only : exit_glc, sigAbort
    !
    ! !ARGUMENTS:
    !
    ! !LOCAL VARIABLES:
    integer(int_kind) :: nml_error   ! namelist i/o error flag

    namelist /glc_override_nml/ enable_frac_overrides, &
         decrease_override_delay, increase_override_delay, rearrange_override_delay, &
         decrease_frac, increase_frac, rearrange_freq
    
    character(len=*), parameter :: subname = 'read_namelist'
    !-----------------------------------------------------------------------
    
    ! Initialize namelist inputs
    enable_frac_overrides = .false.
    decrease_override_delay = 0
    increase_override_delay = 0
    rearrange_override_delay = 0
    decrease_frac = 0._r8
    increase_frac = 0._r8
    rearrange_freq = 0

    ! Read namelist
    if (my_task == master_task) then
       open(nml_in, file=nml_filename, status='old', iostat=nml_error)
       if (nml_error /= 0) then
          nml_error = -1
       else
          nml_error =  1
       endif
       do while (nml_error > 0)
          read(nml_in, nml=glc_override_nml,iostat=nml_error)
       end do
       if (nml_error == 0) close(nml_in)
    end if

    call broadcast_scalar(nml_error, master_task)
    if (nml_error /= 0) then
       call exit_glc(sigAbort,'ERROR reading glc_override_nml')
    end if

    ! Write namelist settings
    if (my_task == master_task) then
       write(stdout,blank_fmt)
       write(stdout,ndelim_fmt)
       write(stdout,blank_fmt)
       write(stdout,*) ' GLC Override Frac:'
       write(stdout,blank_fmt)
       write(stdout,*) ' glc_override_nml namelist settings:'
       write(stdout,blank_fmt)
       write(stdout, glc_override_nml)
    end if

    ! Send namelist settings to all procs
    call broadcast_scalar(enable_frac_overrides, master_task)
    call broadcast_scalar(decrease_override_delay, master_task)
    call broadcast_scalar(increase_override_delay, master_task)
    call broadcast_scalar(rearrange_override_delay, master_task)
    call broadcast_scalar(decrease_frac, master_task)
    call broadcast_scalar(increase_frac, master_task)
    call broadcast_scalar(rearrange_freq, master_task)

  end subroutine read_namelist

  !-----------------------------------------------------------------------
  function frac_overrides_enabled() result(enabled)
    !
    ! !DESCRIPTION:
    ! Return true if glc fraction overrides are enabled in this run
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    logical :: enabled  ! function result
    !
    ! !LOCAL VARIABLES:
    
    character(len=*), parameter :: subname = 'frac_overrides_enabled'
    !-----------------------------------------------------------------------
    
    enabled = enable_frac_overrides

  end function frac_overrides_enabled

  !-----------------------------------------------------------------------
  subroutine do_frac_overrides(ice_covered, topo, ice_sheet_grid_mask)
    !
    ! !DESCRIPTION:
    ! Do all overrides of glc fraction
    !
    ! !USES:
    use shr_log_mod     , only : errMsg => shr_log_errMsg
    !
    ! !ARGUMENTS:
    real(r8), intent(inout) :: ice_covered(:,:)
    real(r8), intent(inout) :: topo(:,:)
    real(r8), intent(in)    :: ice_sheet_grid_mask(:,:)
    !
    ! !LOCAL VARIABLES:
    
    character(len=*), parameter :: subname = 'do_frac_overrides'
    !-----------------------------------------------------------------------

    call apply_increase_frac(ice_covered, topo, ice_sheet_grid_mask)
    call apply_decrease_frac(ice_covered, topo, ice_sheet_grid_mask)
    call apply_rearrange_freq(topo, ice_sheet_grid_mask)

  end subroutine do_frac_overrides

  !-----------------------------------------------------------------------
  subroutine apply_increase_frac(ice_covered, topo, ice_sheet_grid_mask)
    !
    ! !DESCRIPTION:
    ! Apply increase_frac to ice_covered
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    real(r8), intent(inout) :: ice_covered(:,:)
    real(r8), intent(in)    :: topo(:,:)
    real(r8), intent(in)    :: ice_sheet_grid_mask(:,:)
    !
    ! !LOCAL VARIABLES:
    real(r8) :: increase_topo_threshold
    integer(int_kind) :: increase_time_since_baseline

    character(len=*), parameter :: subname = 'apply_increase_frac'
    !-----------------------------------------------------------------------

    increase_time_since_baseline = time_since_baseline(increase_override_delay)
    if (increase_time_since_baseline > 0) then
       ! When increase_time_since_baseline * increase_frac is 0, we'll set elevations >= max_height
       !   to ice_covered = 1
       ! When increase_time_since_baseline * increase_frac is 1, we'll set all elevations >= 0 to
       !   ice_covered = 1
       ! In between those times, we'll use an intermediate threshold
       increase_topo_threshold = (1._r8 - increase_time_since_baseline * increase_frac) * max_height
       increase_topo_threshold = max(increase_topo_threshold, 0._r8)
       increase_topo_threshold = min(increase_topo_threshold, max_height)

       where (ice_sheet_grid_mask > 0._r8 .and. topo >= increase_topo_threshold)
          ice_covered = 1._r8
       end where
    end if

  end subroutine apply_increase_frac

  !-----------------------------------------------------------------------
  subroutine apply_decrease_frac(ice_covered, topo, ice_sheet_grid_mask)
    !
    ! !DESCRIPTION:
    ! Apply decrease_frac to ice_covered
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    real(r8), intent(inout) :: ice_covered(:,:)
    real(r8), intent(in)    :: topo(:,:)
    real(r8), intent(in)    :: ice_sheet_grid_mask(:,:)
    !
    ! !LOCAL VARIABLES:
    real(r8) :: decrease_topo_threshold
    integer(int_kind) :: decrease_time_since_baseline
    
    character(len=*), parameter :: subname = 'apply_decrease_frac'
    !-----------------------------------------------------------------------

    decrease_time_since_baseline = time_since_baseline(decrease_override_delay)
    if (decrease_time_since_baseline > 0) then
       ! When decrease_time_since_baseline * decrease_frac is 0, we'll set elevations < 0 to
       !   ice_covered = 0
       ! When decrease_time_since_baseline * decrease_frac is 1, we'll set all elevations < max_height
       !   to ice_covered = 0
       ! In between those times, we'll use an intermediate threshold
       decrease_topo_threshold = (decrease_time_since_baseline * decrease_frac) * max_height
       decrease_topo_threshold = max(decrease_topo_threshold, 0._r8)
       decrease_topo_threshold = min(decrease_topo_threshold, max_height)

       where (ice_sheet_grid_mask > 0._r8 .and. topo < decrease_topo_threshold)
          ice_covered = 0._r8
       end where
    end if
    
  end subroutine apply_decrease_frac

  !-----------------------------------------------------------------------
  subroutine apply_rearrange_freq(topo, ice_sheet_grid_mask)
    !
    ! !DESCRIPTION:
    ! Apply rearrange_freq to topographic heights.
    !
    ! Example: if rearrange_freq = 3, then for the first 3 days, there will be no
    ! rearrangement, for the next 3 days (days 4-6) the topographic heights will be
    ! rearranged, for the next 3 days (days 7-9) the topographic heights will be back to
    ! normal, etc.
    !
    ! If rearrange_frac is <= 0, no rearrangement is done.
    !
    ! !USES:
    !
    ! !ARGUMENTS:
    real(r8), intent(inout) :: topo(:,:)
    real(r8), intent(in)    :: ice_sheet_grid_mask(:,:)
    !
    ! !LOCAL VARIABLES:
    integer(int_kind) :: num_intervals ! number of intervals of rearrange_freq
    integer(int_kind) :: rearrange_time_since_baseline

    character(len=*), parameter :: subname = 'apply_rearrange_freq'
    !-----------------------------------------------------------------------

    rearrange_time_since_baseline = time_since_baseline(rearrange_override_delay)
    if (rearrange_time_since_baseline > 0 .and. rearrange_freq > 0) then
       ! num_intervals will be 0 in the first interval, 1 in the second, etc.
       num_intervals = rearrange_time_since_baseline / rearrange_freq

       ! Rearrange topographic heights half the time
       if (modulo(num_intervals, 2) == 1) then
          where(ice_sheet_grid_mask > 0._r8)
             topo = max(0._r8, max_height - topo)
          end where
       end if
    end if
    
  end subroutine apply_rearrange_freq


  !-----------------------------------------------------------------------
  integer(int_kind) function time_since_baseline(delay)
    !
    ! !DESCRIPTION:
    ! Return time (days) since the baseline for adjustments (based on delay)
    !
    ! !USES:
    use glc_time_management, only : elapsed_days_init_date
    !
    ! !ARGUMENTS:
    integer, intent(in) :: delay  ! number of days to delay before starting adjustments
    !
    ! !LOCAL VARIABLES:
    
    character(len=*), parameter :: subname = 'time_since_baseline'
    !-----------------------------------------------------------------------
    
    time_since_baseline = elapsed_days_init_date - delay

  end function time_since_baseline




end module glc_override_frac
